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Abstract 

This study seeks a numerical algorithm which optimizes frequency precision for the 
damped sinusoids generated by the nonresonant LITA technique. It compares computed 
frequencies, frequency errors, and fit errors obtained using five primary signal analysis methods. 
Using variations on different algorithms within each primary method, results from 73 fits are 
presented. Best results are obtained using an AutoRegressive method. Compared to previous 
results using Prony’s method, single shot waveform frequencies are reduced ~0.4% and 
frequency errors are reduced by a factor of ~20 at 303K to ~ 0.1%. We explore the advantages of 
high waveform sample rates and potential for measurements in low density gases. 


1. Introduction 

Interest in the nonresonant Laser- Induced Thermal Acoustics (LITA) optical diagnostic 
comes from its capability to make high-accuracy (typically 1%) measurements of temperature[l], 
sound speed[2], velocity[3], and pressure[4] during a single laser pulse. There are few 
instrumental methods combining this level of accuracy and temporal resolution. LITA data can 
be used to advance our understanding in a wide variety of scientific fields. 

Our nonresonant LITA apparatus and the theory behind it have been described previously 
[1]. Briefly, a high-peak-power pulsed (10 ns) laser is split into two beams and crossed at a small 
angle; typically 1 degree. At the focus, the resulting interference pattern creates two counter 
propagating acoustic wave packets by electrostriction. As they counter propagate, they interfere 
constructively and destructively. A second continuous laser beam is Bragg scattered off the 
acoustic wave packets. This creates an intensity modulation in the scattered laser light at a 
frequency (~20 MHz) twice that of the counter propagating acoustic waves (~10 MHz). This 
frequency can be related to sound speed in the gas. Temperature can be calculated if gas 
composition is known. These ultrasonic acoustic waves rapidly decay by acoustic absorption. 
Hence signals reflect the intensity modulation of the counter propagating acoustic waves 
convolved with a damping function created by acoustic absorption. Resulting waveforms 
resemble damped sinusoids. 

2. Previous Analysis Methods 

Data and waveforms analyzed in this study were generated previously (1). Unlike 
reference 1, this study uses raw LITA data not normalized by probe laser temporal variation. 

This weakly perturbs the waveform shape but has negligible effect on computed frequency. This 
approach forces the mathematics to be more robust. 
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In early work on LITA, an analytical expression [4] was derived to describe these 
waveforms. Simplified, it is essentially a sinewave convolved with a decaying exponential which 
must be squared to reproduce the waveform. It is given as 

Y(t) = {A * (exp-k*t) + [B * cos(co* t+ t)i) * (exp- C* t)] } 2 (1) 

where A and B are amplitudes related to the resonant and nonresonant LITA components 
respectively, k is related to gas properties and grating frequency, co is the acoustic frequency of 
interested in this study and related to the sound speed, <)) is a phase factor, and C is related to 
acoustic damping. Initial guesses were required and fits were optimized using a Levenberg- 
Marquardt algorithm. 

Equation 1 is applicable to both resonant and nonresonant LITA. In resonant LITA (i.e. A 
>B), an appropriate molecule such as N0 2 is seeded into the flow and excited to create a strong 
thermal grating and high signal levels on an essentially flat baseline. In our applications with 
large-scale wind tunnel applications, seeding is undesirable. It is expensive, potentially toxic, and 
perturbs the flow. Instead, we pursued nonresonant LITA (i.e. B >A) so the first term in Equation 
1 is presumed negligible. Nonresonant LITA creates acoustic waves in any fluid using 
electrostriction. Unfortunately, nonresonant LITA signals are typically orders-of-magnitude 
weaker than resonant LITA. Signal to noise (S/N) ratios are reduced. Baselines are nonzero and 
often time dependent. Pulse-to-pulse spatial nonuniformities in the laser beam can create a 
waveform where the amplitude of the second peak is larger than the first. All these experimental 
difficulties combine to produce a unique waveform shape that varies on a pulse to pulse basis. 
Hence, nonresonant LITA waveforms only approximately resemble damped sinusoids. 

A variety of gas parameters can be extracted from nonresonant LITA waveforms. Our 
thermometry[l], sound speed[2], velocimetry[3], and pressure [4] applications require only 
waveform frequency. Hence, this study concentrates on minimizing frequency error. Data 
analyzed here was taken at a fixed gas temperature, pressure, and composition. This fixes the 
sound speed in the gas. 

In Reference 1, we recommended Prony’s method. It was developed specifically for 
fitting damped sinusoids. Unlike Equation 1 , it does not require good initial guesses and is 
computationally faster. Results [1] indicate Equation 1 and Prony’s method both produce 
frequencies with ~l-3% single pulse accuracy and precision for nonresonant LITA waveforms. 

In terms of frequency error, these are treated as equivalent methods in this study. 

3. Results and Discussion 

Since Reference 1 was published, new signal analysis methods have become available. 
This study explores several to find the frequency components in the LITA waveforms and 
determines which produce the lowest fit and frequency errors. To eliminate extensive 
programming of complex math routines, we use the AutoSignal (www.systat.com)) software 
package. It is designed specifically to model damped and undamped sinusoids. We select five 
main signal analysis methods (Lourier Transform, AutoRegressive, Prony, Minimum Variance, 
and EigenAnalysis) and test different algorithms associated with each one. The AutoSignal 
manual contains a discussion of advantages of each method, waveform s to which they are 
applicable, and literature references for further reading. Decisions for rejecting some algorithms 
are based on those discussions. AutoSignal abbreviations are adopted here and identified in Table 
1 . 

The errors in this study are based on comparison between experimental and computed 
waveforms. Best results have errors which are sufficiently low that their accuracy cannot be 
verified using current LITA instrumentation. We emphasize that this work is a study in computed 
precision not accuracy. 
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3.1 Waveform Temporal Optimization 

Typical nonresonant LITA waveforms are shown in Figure 1. Initial signal growth is 
dominated by electrostriction. It would be difficult to model this growth even if data were 
available. Acoustic decay of the counter propagating acoustic waves dominates near the first 
peak. After 9-10 cycles, background signal dominates. We expect optimal fit results between 
these locations. Rather than attempt to model these complex effects, we choose the simpler 
solution of selecting the intermediate points which produce minimum computed frequency errors. 




Figure 1. Typical single-laser-pulse damped sinusoidal waveforms produced by nonresonant 
LITA at 303K and 652K. Data is taken from Reference 1. Waveform contains 25 data points per 
cycle at 303K, contains 285 total data points, and is uniformly sampled in 2 ns increments. 
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Figure 2. Computed frequency (MHz, left) and frequency error (%, right) (rms frequency error / 
frequency * 100) as a function of start time for 303 K data in Figure 1 using an AutoRegressive 
spectrum algorithm (Data Svd FB, order 50, signal subspace 4) and non-linear optimization. First 
peak in Figure 1 occurs in Figure 2 at t=0 ns. Second peak occurs at t=50 ns. 

Computed frequency and frequency error results as we sequentially removed initial 
points from each waveform in 4 ns steps are shown in Figure 2. Errors are minimal after initial 
baseline is removed and signals are > 15% of the first peak value. For ease of analysis, we 
discard all data prior to the top of the first peak in Figure 1. We designate it as t=0 in Figure 2 
and fit start in Figure 1. A similar study was performed for the trailing data. They show a very 
weak minimum for both temperatures near the top of the second to last visible peak. However, 
frequencies vary < 0.2% and fit uncertainties < 8% as trailing points are removed. For ease of 
analysis, we discard all data after the second to last visible peak. Resulting fit start and end points 
are shown in Figure 1. Similar result was found for 652K data in Figure 1. 

3.2 Algorithm Comparison 

Algorithm comparison using data in Figure 1 at 303 K is shown in Table 1. Using five 
main signal analysis methods and variations on multiple algorithms within each, we present 
results from 73 fits. Extracting the frequency of damped sinusoids with high precision and 
accuracy is an old and unsolved problem in mathematics. Modem publications recognize the 
AutoRegressive spectrum method (hereafter designated AR) as the best frequency estimator. 
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Table 1. Frequencies, Frequency Errors, and Fit Errors Computed using Single Pulse LITA Data 
in Figure 1 at 303K. Frequency Error (%) = Standard Deviation / Frequency * 100. Dashed Line 
Indicates no Convergence. 

SS = Signal Subspace SVD = Singular Value Decomposition 
NE = Normal Equations 

F, B, FB = Forward, Backward, Forward and Backward Prediction Algorithm 


Fourier Spectrum 


Transform 

Freq, MHz Freq Error, % 

Best Exact N 

20.580 

1.66 


FFT Radix 2 

20.268 

1.60 


Prime Factor 

20.729 

1.68 


Mixed Radix 

20.580 

1.66 


Chirp-Z 

20.580 

1.66 


AutoRegressive Spectrum - All use SS=4, no non-linear optimization 

Algorithm 

Order 

Freq, MHz 

Fit Error, % 

AutoCorr 

40 

20.639 

2.24 

Burg 

30 

20.023 

1.40 

Nrml F 

40 

20.879 

2.15 

Data F 

30 

20.244 

1.50 

Data B 

30 

20.312 

1.78 

Data FB 

30 

20.306 

1.40 

Data SVD F 

50 

20.091 

7.45 

Data SVD B 

50 

20.176 

8.22 

Data SVD FB 

30 

20.218 

3.18 

Data SVD FB 

40 

20.186 

5.35 

Data SVD FB 

50 

20.177 

5.26 optimal freq. fit 

Data SVD FB 

60 

20.180 

7.89 

Data SVD FB 

70 

20.171 

9.39 

Data SVD FB 

80 

20.162 

9.29 

Non-Linear Optimization - Least Squares Minimization Method 

Algorithm 

Order 

Freq, MHz 

Freq. Error, % 

Data SVD FB 

40 



Data SVD FB 

50 

20.169 

0.096 optimal fit 

Data SVD FB 

60 

20.169 

0.096 

Data SVD FB 

70 

20.169 

0.096 

Data SVD FB 

80 
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Initial Algorithm: Data SVD FB, Order 50, SS=4 
Non-Linear Optimization Parameters 

Model: Exp. Damped Sine, 2 frequency components, 
unshared frequency and phase 




Freq. 

Freq. 

Minimization Method 

MHz 

Error, % r 2 

Least Squares 


20.169 

0.096 0.990 optimal fit 

Least Abs Dev 


20.131 

0.11 0.988 

Lorentzian 


20.097 

0.16 0.973 

Pearson 7 


20.153 

0.11 0.985 

Pronv Spectrum - All use SS= 

2 


Algorithm 

Order 

Freq, MHz 

Freq. Error, % 

Damped 

10 

20.417 

2.3 

Damped 

15 

20.252 

2.1 

Damped 

20 

20.210 

2.3 

Damped 

30 

20.245 

2.3 - Ref. 1 result 

Damped 

40 

20.089 

3.4 

Damped 

50 

20.259 

3.4 

Damped SVD 

30 

17.802 

39.5 

Damped SVD 

40 

33.336 

60.7 

Damped SVD 

50 

19.799 

3.3 

Damped SVD 

60 

20.948 

3.4 

Damped SVD NE 

30 

17.964 

35.6 

Damped SVD NE 

40 

17.792 

12.2 

Minimum Variance Spectrum 



Algorithm 

Order 

Freq, MHz 

Freq. Error, % 

SumArBurg 

30 

20.020 

1.6 

SumArBurg 

40 

19.989 

1.6 

SumArDataF 

30 

20.203 

1.6 

SumArDataF 

40 

20.203 

1.6 

SumArDataB 

30 

20.294 

1.6 

SumArDataB 

40 

20.294 

1.6 

SumArDataFB 

30 

20.264 

1.6 

SumArDataFB 

40 

20.264 

1.6 

Musicus 

30 

20.020 

1.6 

Musicus 

40 

19.989 

1.6 

Capon-Kay 

30 

20.203 

1.6 

Capon-Kay 

40 

20.203 

1.6 
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EigenAnalysis Spectrum - All use SS=4 


Algorithm 

Order Freq, MHz 

Freq. Error, % 

MUSIC Fwd 

30 

20.228 

1.6 

MUSIC Fwd 

40 

20.192 

1.6 

MUSIC Fwd 

60 

20.161 

1.6 

MUSIC Fwd 

100 

20.147 

1.6 

EigVec Fwd 

30 

20.229 

1.6 

EigVec Fwd 

40 

20.195 

1.6 

EigVec Fwd 

60 

20.156 

1.6 

EigVec Fwd 

100 

20.066 

1.6 

MUSIC FB 

30 

20.211 

1.6 

MUSIC FB 

40 

20.200 

1.6 

MUSIC FB 

50 

20.179 

1.6 

MUSIC FB 

60 

20.161 

1.6 

MUSIC FB 

80 

20.171 

0.32 

MUSIC FB 

100 

20.152 

0.33 

EigVec FB 

30 

20.218 

1.6 

EigVec FB 

40 

20.200 

1.6 

EigVec FB 

45 

20.186 

0.32 

EigVec FB 

50 

20.179 

0.32 

EigVec FB 

60 

20.152 

0.32 

EigVec FB 

80 

20.173 

0.32 

EigVec FB 

100 

20.156 

0.32 
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In the AR method, if a model can be fit to a data stream, it can be transformed into the 
frequency domain instead of the data upon which it is based. Hence, it produces a continuous and 
smooth spectrum. In an AR model, a value at time t is based on a linear combination of prior 
values (Forward prediction), combination of subsequent values (Backward prediction) or both 
(Forward-Backward prediction). This is a linear model which gives rise to rapid and robust 
computations. The singular value decomposition (SVD) improves results using in-situ noise 
filtration. 

AR has several algorithm variations including Data, AutoCorr, Burg, and Normal (Nrml). 
Most can be run using F, B or FB prediction. Each of these can vary order creating a large 
number of combinations. Discussions of these methods are provided to minimize the number of 
variations tested. 

AutoCorr (sometimes referred to as Yule-Walker method) is common in most statistical 
packages but has the least frequency resolution. The Burg Algorithm (Max Entropy Method) is 
most widely known but designed for harmonic detection. The NRML algorithm uses a least- 
squares normal equations approach. Sums incur loss of precision which can lead to reduced 
frequency resolution. The DATA algorithm processes the full F, B, or FB prediction data matrix. 
This approach has no loss of precision but requires extensive computation. It is widely 
recognized as the most accurate of AR processing procedures. The DATA FB algorithm is also 
known as the modified covariance method. The DATA SVD FB algorithm is known as the 
Principal Component AutoRegressive method (PCAR). The SVD procedure uses an 
eigendecomposition to discard the noise eigenmodes when the AR coefficients are computed. The 
combination of DATA, SVD, and FB produces an algorithm with extraordinary frequency 
estimation. 

AutoSignal software has a feature called AR spectrum with algorithm comparison. Using 
data in Figure 1, it predicts best results using the Data SVD FB algorithm over all others 
assuming two frequencies in the dataset (signal subspace = 4). Signal subspace is twice the 
number of spectral frequency components. One is for the real and one for the imaginary 
component. Optimal model order determination is an art at best. As model order increases, more 
of the data trends are incorporated into the model reducing uncertainty. If an order is selected 
which is greater than the optimum order, this adds noise which reduces uncertainty. AutoSignals 
order exploration feature predicts an optimum order between 40 and 50. This corresponds to two 
cycles in Figure 1. However, we still chose to explore errors as a function of model order before 
selecting an optimum order. We fit the spectrum over its full range with an adaptive n. We set 
added noise to its minimum value of 300 dB. This adds a fractional noise of IE- 15. AR results 
are listed in Table 1. 

There are two types of errors listed in Table 1. Most algorithms report a frequency and 
frequency error. AR reports fit errors. These are rms errors between the data and the computed fit 
results. These are the exception in Table 1. Nonlinear optimization fits report frequency errors. 
Care must be exercised when comparing these results. Particularly with AR, a high fit error does 
not imply a high frequency error. In fact, we show that AR produces the best frequency error 
results but with poor fit errors. 

All results, particularly fit errors, from any algorithm can be further refined using 
nonlinear optimization (NLO). Autosignal uses an exponentially damped sinusoid model. 

Also, our nonresonant LITA waveforms have a time dependent baseline which we find can be 
approximated by a weakly damped sine wave. So, instead of using Equation 1, we use Equation 
2 below. 


Y(t) = Ai * sin(ooi* t+ <|)i) * (exp- Bi* t) + A 2 * sin(eo 2 * t+ <j> 2 ) * (exp- B 2 * t) (2) 



It is a linear combination of two exponentially damped sinusoids. The first damped sinusoid 
models a low frequency oscillation which is a combination of the LITA baseline and the temporal 
profile of our probe laser. The second damped sinusoid models a high frequency oscillation 
representing the acoustic wavepackets. Both terms have independent phases and damping. 
Outputs from the AR method (DATA SVD FB, Order 50, SS=4) are used as initial guesses. 

LITA ordinate data typically spans ~1 order-of-magnitude. We expect linear-least-squares 
minimization to provide best results. Results in Table 1 confirm this and demonstrate no need for 
any of the three robust minimization approaches. We studied output with NLO as a function of 
order. We define optimal order as a compromise between minimal errors and minimal 
computation time. Hence, we select the lowest order possible. Table 1 shows optimal results 
using an order of 50. Since NLO should provide nearly optimum results, we conclude Figure 1 
frequency is 20.169 MHz with 0.096% frequency error. Thi s is used as the benchmark for 
comparison of frequency and error results between the different algorithms. 

We distinguish between two types of results: those which come from a specific algorithm 
such as AutoRegressive and those which come from the nonlinear optimization routine. Virtually 
all algorithm results can be selected as input guesses to the NLO routine. If functioning properly, 
NLO results will converge to identical values independent of input guesses. However, different 
algorithms or improper order selection can provide poor guesses producing nonconvergence. 

Optimal algorithms in Table 1 are not necessarily those with the lowest errors. Rather, it 
is those whose frequencies change least after nonlinear optimization and converge using a 
minimum number of iterations. Hence, they provide the best guesses to the optimization methods 
and are most likely to produce convergence. Within the AR method. Table 1 results indicate Data 
SVD FB is the optimal algorithm choice. This conclusion is supported by the math behind the 
method, discussions in the AutoSignal manual, and several trial fits with a few other algorithms 
(Burg, Autocorr, Nrml) using a few parameter variations (F, B, and FB) not included in Table 1. 
Using order = 50 and signal subspace = 4 (i.e. 2 frequency components), it computes 20.177 MHz 
with 5.26% fit error). Nonlinear optimization only refines this to 20.169 MHz with 3.06% fit 
error and 0.096% frequency error. There is a 0.0397% difference in the computed frequencies. 

Autoregression computes a fit error but not a frequency error. The nonlinear 
optimization routine computes both. As shown above, both methods produce essentially identical 
frequency results. We conclude the error computed using nonlinear optimization can be assigned 
to the frequency computed using Autoregression with the DATA SVD FB algorithm with order 
50 and signal subspace 4. 

Fourier transforms are standard methods for undamped sinusoids. However, they have 
limited frequency resolution since computed frequencies are equally spaced and fixed in number. 
They are known to be inaccurate when applied to damped sinusoids. Table 1 results using five 
nearly equivalent Fourier algorithms show frequencies are 0.5-2. 8% higher than NLO results with 
~17x larger frequency errors. 

Prony’s method fits a linear combination of complex damped sines (AutoSignal 
approach) or complex damped exponential decays (ref. 1 approach) to uniformly sampled data. 
Damped algorithms are used since undamped methods are inappropriate for these waveforms. 
Results are presented as a function of order using SVD and NE algorithms. We were surprised 
that using SVD did not improve the frequency errors. SVD NE is designed for rapid computation 
of large datasets. It is known to produce results inferior to SVD and is very sensitive to noise. 

Fits are weakly sensitive to order (defined for Prony’s as number of exponential decays fitted). 
Best results in Table 1 (Damped, order=30 SS=2) show frequencies (20.245 MHz) are ~ 0.38% 
higher than NLO results with 24x larger frequency errors. Since results are based on signal 
subspace of 2, they cannot be used as initial guesses for the NLO routine. Results from Prony’s 
method used in reference 1 produce a frequency of 20.237 MHz and frequency error of 2.3%. 
These are nearly identical to Table 1 results. 
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Minimum Variance(MV) and EigenAnalysis(EA) methods are spectral estimators 
designed to discover harmonics in a waveform and graphically compare component powers. 

Their frequency resolution is typically slightly better than FFT’s. Best MV results in Table 1 
show frequencies within 0.17% of NLO results. Flowever, most frequencies are -0.6-0. 8% 
higher than NLO results with 17x larger frequency errors. MV has no provision for order 
variation. NLO will not converge with these inputs. 

EA results in Table 1 show frequencies nearly identical with NLO results. Music Fwd 
and EigVec Fwd have 16x larger frequency errors and wont converge. Music FB and EigVec FB 
have only 3x larger frequency errors and NLO converges. Fligher orders may improve results. 
Autosignal EA routines are limited to order 100. 

Based on Table 1 results, we conclude NLO using Equation 2 produces the most precise 
fit errors, frequency, and frequency errors (-0.1%). This requires initial guesses supplied by the 
AR spectral method using Data SVD FB algorithm with order = 50 and SS=4. If only frequency 
information is required, the AR method using Data SVD FB algorithm provides results nearly 
equivalent to the NLO results but with minimal computations and no initial guess requirement. 
NLO requires - 0.1 sec using a 3.2 GFlz Pentium 4. EA is also viable but with increased 
frequency uncertainty(~0.3%). All other algorithms including Prony’s produce inferior 
results(~1.5-3% frequency errors). 

3.3 Mathematics of Selected Algorithms 

3.3.1 Prony’s Method Spectral Analysis Algorithms 

The waveform to be analyzed is Y(t) and consists of N evenly spaced samples. For 
analysis using Prony’s method, it can be expressed either as the sum of complex conjugate and 
real exponentials 

m 

Y(t)= ^a(j) * exp(-z co(j)t + ®(y)) * exip(-k(j)t ) (3a) 

j = i 

or as the sum of exponentially damped sinusoids 

m 

Y(t)= ^a(j) * exp(-k(j)t)* cos(co(j)t + ®(y)) (3b) 

j = i 

Flere, m is the order and refers to the number of equations which will be summed to 
describe the waveform. For a given equation, a(j) is the amplitude, k(j) is the damping 
coefficient, co(j) is the frequency in hertz, and i|)(j) is the phase in radians. Prony’s method 
computes the complex parameters of each equation and then provides the four real 
parameters(a(j), k(j), eo(j), and ())(])) for each equation. 

We fit Y(t) data in Reference 1 to Equation 3a and use m= 40. AutoSignal uses Equation 
3b. AutoSignal fits in Table 1 are based on one frequency component (signal subspace of 2). 
Both approaches produce equivalent results. 
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3.3.2 AutoRegressive Spectral Analysis Algorithm 


The AR model is defined as follows 

Y(k)= f j a(j)*x(k + j) , k = p...l (4) 

j = i 

Y(k)= ^a(y) *x(k-j) , k = p+l, ...N 

j = i 

Where x is a data series of length N and a is AutoRegressive parameter array of order p. 
AutoSignal uses a positive sign convention(linear prediction) for the AR coefficients. This model 
is defined as backward prediction for the first p values and forward prediction for the remaining 
N-p values. 

3.4 Comparison Between Autoregressive (AR), Autoregressive with Non Linear 
Optimization (AR+NLO), and Prony’s methods 

3.4.1 Single Waveform Results 

Using data in Figure 1, we examine both graphically and numerically the results from 
AR, AR with NLO, and Prony’s method. We begin by comparing power spectra results between 
AR and Prony’s method. Results for 303K data in Figure 1 are shown in Figure 3. AR produces 
both a sharper frequency distribution and better S/N (i.e. 45 dB versus at best 30 dB) for Prony’s. 
Clearly, AR can generate a more precise frequency and frequency error. 




Figure 3. Power spectrum in dB using Prony’s method (left) (Damped, order=30, SS=2) and 
AutoRegressive method (right) (Data SVD FB order=50 SS=4) for 303K data in Figure 1. 

The results of the AR method provide more accurate inputs for the NLO routine. AR 
with NLO computed waveforms using Equation 2 are shown (Figure 4A). These are summed and 
overlaid with Figure 1 data (Figure 4B). Fit results using AR, AR with NLO, and Prony’s 
method are given in Table 2. Component 1 is defined as the low frequency baseline. Component 
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2 is defined as the high frequency acoustic wavepackets. Visually, the frequency accurately 
matches the data throughout the entire waveform. This is in contrast to Figure 3 in Reference 1 
using Prony’s method. It appears to accurately match the frequency early in the waveform but 
predicts a slightly higher frequency with increasing time for both the 292 and 652K data. Results 
in Table 2 using Prony’s method produce a frequency (20.245 MHz) which is 0.38% higher than 
AR with NLO(20. 169 MHz) and AR without NLO. Frequency errors comparing Prony’s and AR 
with NLO in Table 1 differ by a factor of ~18. 

Comparing AR with AR with NLO, Table 2 shows nearly identical frequencies, phases, 
and powers. However, AR results for amplitude are inferior to those of NLO. This produces 
larger AR fit errors. We conclude AR with NLO produce the optimal fit results. If only 
frequency information is required, we recommend the AR method. Since no NLO processing is 
required, it is computationally faster. It requires no initial guesses and produces equivalent 
frequency results. We note that the AR spectrum is produced by modeling the waveform. Hence 
it is not artificially constrained as is Equation 2. Therefore, one could argue that it may produce 
frequency results which are more accurate than the NLO routine. 

Our current LITA apparatus has insufficient accuracy to experimentally validate these 
numerical results. Since both AR and AR with NLO use different math methods to compute their 
results, this provides confidence in the frequencies and their associated errors. Using this fact, 
Prony’s method frequencies and errors, and the Prony’s graphical fits in reference 1, we conclude 
Prony’s method is producing the poorest fit results. 
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Table 2. Comparison Between AutoRegressive (AR), AutoRegressive with Non Linear 
Optimization (AR+NLO) Fit Parameters (Equation 2), and Prony’s Methods Numerical Results 
using Waveform in Figure 1 at 303K. AR + NLO Results are Used to Generate Waveform s in 
Figure 4B. 





Freq, 


Damping 

% 

Fit 

Method 

Sinusoid Amplitude MHz 

Phase 

MHz 

Power Error 

AR±NLO 

Component 1 

205 

0.479 

0.387 

3.97 

91.3 



Component 2 

86 

20.169 

0.46 

5.54 

8.7 

3.06% 

AR 

Component 1 

109 

0.277 

2.04 

— 

94.4 



Component 2 

22 

20.177 

1.43 

— 

5.6 

5.26% 

Prony 

Component 1 

64 

20.245 

1.46 

4.91 

100. 

— 


3.4.2 Results Using 100 Waveforms 

To insure results in Tables 1 and 2 are not limited only to Figure 1 data, we analyzed 100 
waveforms from reference 1 using AR, AR with NLO, and Prony’s methods at both 303K and 
652K. Single shot results of frequency, frequency error and their associated histograms are 
presented in Figures 5-7. Figure 1 data is the first shot in this series. Average results and 
uncertainties are presented in Table 3. 

In Figure 5, frequencies at 303K as a function of laser shot tend to be symmetrically 
distributed about 20.3 MHz. Figure 5 histogram shows shots are tightly grouped within ±0.3% 
(57 shots at 20.30 ±0.06 MHz). Figure 6 histogram shows errors are tightly grouped (83 shots 
between 0.08% and 0.10% error). There are 17 shots with errors > 0.1%. We inspected 
individual shots but could find no correlation between frequency and frequency error. 

Frequencies and frequency errors using Prony’s method are shown in Figure 7 and Table 
3. Average frequency is 20.366 MHz with 1.76% error. 

Using all data in Figures 5 and 7, we computed an unweighted mean and standard 
deviation. A similar analysis was performed with the frequency errors in Figure 6 and Figure 7. 
Results are shown in Table 3 for 303 and 652K data. Prony’s frequencies are centered near 
20.366 MHz with average error of 1.76% while AR with NLO frequencies are centered at 20.305 
MHz with 0.095% error. These algorithms produce average frequencies which differ by ~0.4% 
and average frequency errors which differ by a factor of ~18. These results are similar to those 
above for the single waveform. They confirm our conclusions and demonstrate that the AR 
method is more precise when computing the frequency and frequency error over this ensemble of 
spectra than Prony’s method. However, the unweighted standard deviation of the mean of 100 
frequencies produced using AR analysis did not improve relative to Prony’s method. We 
attribute this to pulse-to-pulse frequency fluctuations in our apparatus. 

Frequency fluctuations can arise from temperature fluctuations in the oven and 
fluctuations in the LITA apparatus. Results in Figure 5 are uniformly distributed about 20.3 MHz 
as a function of time. This suggests oven temperature was essentially constant during the 1 0 
seconds required to acquire 100 shots. This is supported by the following facts. We used a 
sealed cell isolated from room air currents. Cell is surrounded by oven weighing 165 lbs (i.e. a 
large thermal mass). For data at 303K, no power was applied to oven. The result is a slow cool 
down to room temperature and a presumed constant temperature during a 10 second data 
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acquisition time. We conclude the main contribution to frequency fluctuations in Figure 5 are 
fluctuations in the LITA part of our apparatus. 




A? Af v <j? 

Frequency, MHz 


Figure 5. Autoregressive method with nonlinear optimization (AR with NLO) computed 
frequencies (left) and histogram (right) using 303K data for 100 laser shots. 57 shots are at 20.30 
± 0.06 MHz. 



Figure 6. Autoregressive method with nonlinear optimization (AR with NLO) computed 
frequency errors in % (left) and histogram (right) using 303K data for 100 laser shots. 83 shots 
have between 0.08 and 0.10% error. 
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Figure 7. Prony’s method computed frequencies(left) and frequency errors (%, right) using 303K 
data for 100 laser shots. (Damped, order=30, SS=2) 

We also computed the average fit frequency and frequency error using a weighted 
average. Results are shown in Table 3 for 303 and 652K. Weighting reduces the frequency 
uncertainty of both methods by an order of magnitude. Essentially, this is a simple method to 
discard outlier results. By comparison, results in Reference 1 used residuals, noise, along with 
power spectral peak width and area as criteria to reject outliers. Typically 90 shots passed. Here, 
based on the histogram in Figure 6, 17 shots with errors above 0.1% are weakly weighted. Still, 
with AR’s higher precision, 83 shots contribute to the weighed average. AR with NLO results 
indicate an average weighted frequency error of 0.009%. 

To further study differences in the frequencies obtained using these algorithms, we 
calculated for 100 shots on a shot by shot basis, individual percent frequency differences 
between Prony’s and AR with NLO (i.e. (P-NLO)/NLO *100) and frequency differences between 
AR with NLO and AR without NLO (i.e. (AR-NLO)/NLO *100). Results and histograms are 
given in Figures 8 and 9. On a shot by shot basis, Prony’s frequencies are 0.3-0. 6% higher versus 
AR with NLO while AR without NLO are < 0.10% of the AR with NLO results. 

The results of this study indicate thermometry results in Reference 1 suffer minor errors 
due to Prony’s method inadequacies. However, errors are well within the stated uncertainties of 
that work. Also, results are based on sound-speed ratios referenced to room temperature results. 
This tends of minimize small frequency errors. 
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Table 3. Frequencies (mean ± la) and Frequency Errors (% error of fits), (mean ± la) 
Associated with Weighted and Unweighted Analysis of 100 LITA Waveform s using AR with 
NLO and Prony Methods. 

% Error 



AR with NLO 

Prony 


Improvement 





Factor 

Unweighted 




T(K) 

Freq, MHz % Error 

Freq, MHz 

% Error 


303 

20.305 ± 0.15 0.095 ±0.027 

20.366 ±0.091 

1.76 ± 0.14 

19x 

652 

30.43 ±0.17 0.178 ±0.042 

30.539± 0.169 

3.15 ± 0.56 

18x 

Weighted 




T(K) 

Freq, MHz % Error 

Freq, MHz 

% Error 


303 

20.29720 ± 0.00183 0.0090 

20.364 ±0.036 

0.177 

20x 

652 

30.4381 ± 0.0051 0.0168 

30.5396 ±0.093 

0.305 

18x 


17 





Frequency Error, % 


Figure 8. Data analysis comparison between Prony’s method(P) and Autoregressive with 
nonlinear optimization (NLO) i.e. (P-NLO)/NLO *100 using computed frequencies(left) and 
histogram(right) for 303K data for 100 laser shots. 




Frequency, MHz 


Figure 9. Data analysis comparison between AutoRegressive(AR) method and AutoRegressive 
method with nonlinear optimization(AR with NLO). (AR-NLO)/NLO* 1 00 using computed 
frequencies(left) and histogram(right) for 303K data for 100 laser shots. 


18 


3.4.3 Discussion 


There are several reasons why AR produces frequency and frequency error results 
superior to Prony’s method. 

First, the AR spectrum is produced by modeling as opposed to other methods such as 
Prony’s which are transforms. It is a linear non-iterative matrix method. The DATA and FB 
options process the lull forward and backward prediction data matrix to produce high accuracy 
results. It is known to have less bias in frequency estimate of spectral components and reduced 
variance in frequency variance over an ensemble of spectra(6). It is particularly well suited for 
the limited number of data points typical of LITA waveforms. Adding the SVD method forces 
robust predictions which are based on signal components, minimally influenced by noise, and 
exceedingly sharp. The adaptive n option uses a Runge-Kutta procedure to integrate the spectrum 
adaptively. This produces a frequency set concentrated near the peaks and an accurate area under 
the spectrum. 

Second, AR results are further refined albeit modestly using nonlinear least-squares 
optimization with a standard Levenberg-Marquardt routine using Equation 2. Fit errors are 
reduced by a factor of 2 while frequencies changes by < 0.1%. Although not derived rigorously 
as Equation 1 , Equation 2 describes the nonresonant LITA waveform produced by our 
experimental apparatus. It is limited to two frequency components and is an 8 parameter fit. 
Prony’s requires 40 equations. This results in less interaction between fit parameters and more 
certain results. 

At first glance, these exceptional results seem unjustified since Equation 2 has a different 
functional form compared to Equation 1. As discussed in section 3.2, the first damped sinusoid 
well describes the combined temporal profiles of our probe laser and LITA baseline. Our probe 
laser is a pulse stretched Alexandrite whose temporal profile compares well with computed 
baseline waveform in figure 4A. By comparison with equation 1, the second damped sinusoid in 
equation 2 should be squared. However, in nonresonant LITA, pulse-to-pulse spatial 
nonuniformities in the pump laser create nonuniform intensities in the counterpropagating 
acoustic wavepackets. Therefore, neither equation 1 nor the second sinusoid in equation 2 are the 
true functional form. The true function is between these extremes and pulse dependent. Hence, 
we consider the use of a damped sinusoid to be a reasonable approximation. Combined, these 
two effects are reasonably well described as the sum of two damped sinusoids. We conclude 
Equation 2 is justified. 

The increased frequency resolution possible using the AR method has revealed a 
curiosity in this LITA data. Two frequency component fit residuals show a clear oscillatory 
behavior. This generally indicates another frequency component exists in the data. To test this 
we refit the data using 3 components (i.e. signal subspace 6). The fits are in a word perfect with 
r2 = 0.999. Data in Figure 1 at 303K produces 19.090 MHz (±0.50%) and 20.728 MHz (± 
0.29%). Amplitude, phase, damping, and power are similar for both frequencies. However, we 
see no physical reason for the existence of a third frequency. Currently, we consider it an artifact 
of the model, the fitting routine, our LITA apparatus, or the low data sampling rate. However, we 
observe nearly symmetrical frequency splitting and similar characteristics of amplitude, phase, 
damping and power. One possibility is that we are seeing the beat frequency between the low- 
frequency oscillation (0.479 MHz) due to the combined baseline and probe laser profile and the 
high frequency oscillation due to the acoustic wavepackets(20.169 MHz). This would produce 
20.169 - 0.479 = 19.69 MHz and 20.169 + 0.479 = 20.648 MHz i.e. nearly our observations. A 
more remote possibility is that we are somehow detecting frequency differences in the counter 
propagating acoustic wavepackets. It is hoped an improved LITA apparatus with a higher 
sampling rate (discussed below) will shed some light on this curiosity. No evidence for a fourth 
frequency component was detected. 
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3.5 Sample Rate Study 

The AR power spectral density function is defined as follows[7] 
a 2 At 

Par (®) = (5) 

p 

{1+ '^ j a(j)*exp(-icokAt)} 2 

j = i 

Here, a is the white noise variance and A t is the sampling interval. This equation shows that fit 
uncertainty improves as sampling interval decreases. Although derived for the Autoregessive 
method, the general concept is applicable to all methods. Hence we undertook a sampling rate 
study to understand how this can further improve fit precision and accuracy. 

Table 4 provides data on frequency error improvement with increased sampling rate 
using the AR analysis method with nonlinear optimization. For this study, we require a highly- 
stable highly-ac curate frequency source. We choose a function generator which produced a 100 
mV undamped sine wave at 20.000 000 MHz with ±5 ppm stability(0.0005% frequency error). 
We acquired 500 ns of data at sampling rates of 0.1, 0.5, 2.5, and 20 GHz using a digital 
oscilloscope (2.5 GHz analog bandwidth, 40 GHz sampling rate). Results in Table 4 show that as 

sample rate is increased by a factor of N, frequency error is reduced by y[N . This is essentially 
Felgetf s advantage. However, it was achieved with sampling rates well above the analog 
bandwidth of the oscilloscope. 

Results in Table 4 at 20 GHz show frequency error (0.0016%) is only a factor of ~3 
greater than the manufacturer’s instrument uncertainty(0.0005%). Hence, the AutoRegressive 
method is capable of this level of accuracy. Note that LITA result errors, albeit for damped 
sinusoids, are nearly two orders-of-magnitude larger. This study provides confidence in the 
precision and accuracy of the fit errors computed for damped sinusoids. 

Table 4. Sample Rate Study using AutoRegressive Method (Data SVD FB, Order 40, Signal 
Subspace 2) with Non-Linear Optimization. 


Sample 

Rate, 

GHz 

Data 

Points 

Computed 

Frequency, 

MHz 

Frequency 
% Error 

Improvement Factors 
a/n Theory Measured 

0.1 

50 

20.000 

0.0191 



0.5 

250 

19.997 

0.0088 

2.2 

2.2 

2.5 

1250 

19.999 

0.0043 

5.0 

4.4 

20.0 

10000 

20.000 

0.0016 

14. 

12. 


3.6 Cycle Reduction Study 

We are interested in applying LITA to supersonic and hypersonic wind tunnels along 
with combustors. These are low-density applications where acoustic damping is large. The result 
is a limited number of cycles and larger fit and frequency errors. Hence we studied errors as a 
function of number of cycles for damped and undamped sinusoids. We used AR with NLO and 
Prony’s method to quantify errors. 
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3.6.1 Damped Sinusoid 


Results of a cycle reduction study comparing AutoRegressive method with nonlinear 
optimization to Prony’s method using Figure 1 data at 303K are shown in Table 5. For AR with 
NLO, they show 0.14% frequency errors with 2 cycles. Data from 1.5 cycles produced a 
frequency but NLO results would not converge. By comparison, Prony’s method produces 
comparable frequency results but with errors which increase from 2% with 9 cycles to 12% with 
1.5 cycles. Prony’s method errors are 20 times larger than AR with NLO over 9 cycles and 57 
times larger than AR with NLO over 2 cycles. We conclude AR or AR with NLO is the superior 
approach for a limited number of cycles. We conclude by using AR or AR with NLO, LITA can 
be a viable diagnostic in low-density applications. 

Table 5. Frequencies and Frequency Errors Computed as a Function of Cycles for Data in Figure 
1 using AutoRegressive Method (DATA SVD FB, SS=4) with Nonlinear Optimization and 
Prony’s method. 



AR with NLO 


Prony’s 



Freq, 

Freq. 


Freq, 

Freq, 

Cycles 

MHz 

% Error 

Order 

MHz 

% Error 

9 

20.168 

0.0983 

50 

20.303 

2.05 

8 

20.170 

0.101 

50 

20.305 

2.27 

7 

20.144 

0.109 

50 

20.303 

2.55 

6 

20.187 

0.106 

80 

20.292 

2.94 

5 

20.179 

0.116 

80 

20.308 

3.46 

4 

20.300 

0.142 

50 

20.321 

4.14 

3 

20.224 

0.138 

50 

20.335 

5.25 

2 

20.295 

0.140 

35 

20.257 

8.10 

1.5 

20.496 


— 

20.403 

12.0 


3.6.2 Undamped Sinusoid 

To further evaluate Table 5 results, we require knowledge of how accurately AR with 
NLO can extract the frequency and error of a known waveform as a function of number of cycles. 
To perform this study, data was simulated using the function generator described above to create 
a 20.000 000 MHz undamped sine wave with 200 mV peak value which fills the digitizer. 
Waveform was sampled using the digital oscilloscope described previously. Results are shown in 
Table 6 for 0.5 and 20 GHz sampling rates 

For 10 cycles at 20 GHz, per cent frequency errors are only 3 times the stated function 
generator uncertainty. Errors for 9 cycle damped sine waves in Table 5 are 52 times higher. 
Although we are comparing damped and undamped sine waves, this study provides confidence in 
the computed frequencies and uncertainties produced using AR with NLO. 

Using 0.5 GHz sampling, fits were stable with three cycles. For two cycles, results were 
unstable (i.e. they seldom converged). Data from one cycle would not converge. For 20 GHz 
sampling, fits were stable to two cycles. Minor instability was encountered with one cycle. 
Repeated attempts would occasionally produce convergence at 0.5 cycles albeit with higher 
uncertainties. Comparing 20 versus 0.5 GHz sampling rate for any specified cycle, we observe 
that % error follows a simple linear dependence. Within a given sampling rate, % error increases 
nonlinearly as number of cycles is reduced. Using an oversimplified model and plotting the 
inverse % error (1/y) versus number of cycles (x), we find at 0.5 GHz and 20 GHz respectively. 
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1/y = (-19.0±2.6) + ( 12.8±0.4)*x 
1/y = (-49.3±6.2) + (52.4±2.8)*x 


( 6 ) 

( 7 ) 


We conclude that using 20 GHz sampling, AR with NLO can determine frequencies with 
errors < 0.05% over 1-2 cycles for undamped sine waves. 

Table 6. Frequencies and Frequency Errors as a Function of Cycles for 20 and 0.5 GHz Sampling 
of a 20.000 000 MHz Simple Undamped Sine Wave using AutoRegressive Spectrum (Data SVD 
FB, Order 50, Signal Subspace 2) with Non-linear Optimization. 



20 GHz sampling 

0.5 GHz sampling 


Freq, Error, 

Freq, 

Error, 

Cycles 

MHz % 

MHz 

% 

10 

20.000 0.0016 

19.997 

0.009 

9 

20.000 0.0019 

19.998 

0.010 

8 

20.000 0.0023 

19.999 

0.012 

7 

20.000 0.0027 

20.000 

0.014 

6 

20.000 0.0034 

19.996 

0.019 

5 

20.001 0.0046 

19.995 

0.023 

4 

20.002 0.0065 

19.997 

0.035 

3 

19.999 0.0097 

20.002 

0.061 

2 

19.991 0.0179 

19.987 

0.093 

1.5 

20.000 0.0275 

20.0 

0.221 

1 

19.989 0.0502 



0.5 

20.3 0.110 



3.7. Signal Reduction Study 




LITA is a nonlinear optical method. Signal levels can vary significantly on a pulse to 
pulse basis. It is useful to study how errors vary as a function of signal level over the digitizer 
range for high sample rates. Unfortunately, our current LITA apparatus is not configured for 20 
GHz sampling. We choose to study this effect using undamped sine waves created using the 
function generator described above. Results at five signal levels spanning 90% down to 4.5% of 
digitizer range are shown in Table 7. Errors are independent of signal level to approximately 
18% of the digitizer range. Below this, they increase by a factor of 2.5 down to 4.5% of the 
digitizer range. By comparison, when damped sinusoid signals similar to those in Figure 1 
dropped below 15% digitizer range, they often could not be fit with Prony’s method. These 
results suggest AR with NLO is better than Prony’s at fitting waveforms with low S/N. 

Table 7. Signal Reduction Study using AutoRegressive Spectrum (Data SVD FB, Order 50, 
Signal Subspace 2) with Nonlinear Optimization of a 20.000 000 MHz Undamped Sine Wave 
using 10 cycles Sampled at 20 GHz. 


Signal, 

Freq, 

Error, 

mV 

MHz 

% 

200 

20.000 

0.0020 

100 

19.999 

0.0015 

40 

19.999 

0.0019 

20 

19.998 

0.0028 

10 

19.999 

0.0048 
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4. Conclusions 


This study demonstrates that the AutoRegressive spectral analysis algorithm with 
nonlinear optimization(AR with NLO) produces optimal frequency results, frequency errors, and 
fit errors for analysis of nonresonant L1TA waveforms with homodyne detection. The nonlinear 
optimization model is a linear combination of two damped sinusoids. Using 0.5 GHz sampling 
rate, this method provides 0.096% single laser pulse frequency error at 303K. This is a factor of 
~20 improvement compared to previous results using Prony’s method. 

If only frequency information is required, this study recommends the AutoRegressive 
spectral analysis method using the Data SVD FB algorithm. For conditions of this study, results 
are optimized using a model order of 40 and signal subspace of 4. Frequencies are equivalent to 
those of the AR with NLO algorithm, require no initial guesses, and require less computation. 
Since the AR spectrum is produced by modeling, it should produce frequency results which are 
equivalent to the NLO routine. Since AR and AR with NLO compute equivalent frequency 
results and frequency errors using different mathematical approaches, they independently verify 
the computed frequencies and frequency precision. 

Results are optimized for parameters specific to these datasets. If any parameter 
particularly sampling rate is altered, relevant sections of this study should be repeated to 
determine optimal fit conditions. Model order can be particularly sensitive to this effect. 

This study demonstrates three advantages of higher sampling rates with undamped 
sinusoids. As sample rate is increased by a factor of N, frequency error is reduced by a/n. Fits 
converge with low errors with a limited number of waveform cycles. Fits converge and errors are 
reduced with low signal levels. Similar results are expected for damped sinusoids. Therefore, 
this study recommends the highest possible sampling rate to reduce errors. 

Results obtained by combining AR analysis with increased sampling rate suggest 
nonresonant L1TA can be applied with high accuracy and precision in low-density environments 
such as combustors and hypersonic wind tunnels. Results suggest single shot errors < 1% are 
possible even when number of waveform cycles are limited to 2 cycles by acoustic damping. 

This study suggests 0.001% or 1 part in 10 5 sound speed measurements are possible by 
combining the AR algorithm for data analysis (0.1% or better single shot frequency and 
frequency error precision), a high-sample-rate data acquisition system ( lOx error reduction using 
50 GHz sampling), and waveform averaging (lOx error reduction using 100 waveform average). 
Therefore, we advocate a long-term development effort with two goals. First, attempt to develop 
a high-stability high-sample-rate laboratory L1TA instrument where most shots have 0.01% or 
better frequency accuracy and precision. Second, use existing methods(8) to build an oven with 1 
mK temperature resolution and thermometer capable of verifying the 0.001% error level. 

LITA measurements are the result of sound propagation in the free field. This is in 
contrast to standard methods such as resonators (spherical and cylindrical) and ultrasonic 
interferometers which are encumbered by boundaries. Uncertainties of the latter are typically 
0.01% or less (9). Results of this study suggest that with additional development, LITA can rival 
these methods for fundamental acoustic measurements without the boundary issues. 
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